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1  Forward 


Our  research  program  has  emphasized  innovative  computations  and  theory.  The  computations  and  theorv 
support  and  enhance  each  other.  We  have  a  coherent  approach  which  depends  upon  abstracting  impor¬ 
tant  mathematical  concepts  and  computational  methods  from  individual  applications  to  a  wide  range  of 
applications  involving  complex  continua.  including  wave  refractions,  flows  in  elastic  and  plastic  media,  and 
complex  fluid  mixing.  We  have  also  developed  adaptive  computational  methods  for  flows  with  discontinu¬ 
ities  and  implemented  these  methods  on  modern  parallel  computers. 

2  List  of  Appendixes 

There  are  no  appendixes  in  this  report. 

3  Report  of  Research 

3.1  Statement  of  Problem  Studied 

3.2  Summary  of  Most  Important  Work 

3.2.1  Discontinuities  and  Adaptive  Computation 

The  front  tracking  method  [9],  [  20]  is  a  computational  method  that  incorporates  explicit  degrees  of  freedom 
to  represent  dynamical  interfaces  and  wave  interaction  fronts.  The  method  is  highly  developed  in  two  space 
dimensions,  and  allows  the  resolution  of  complex,  chaotic  interfaces  between  two  interpenetrating  fluids 
[18].  It  also  allows  for  the  refraction  of  shock  waves  by  interfaces  in  a  number  of  cases. 

Interface  methods  (such  as  front  tracking)  are  the  only  computational  methods  to  duplicate  the  exper¬ 
imentally  correct  growth  rate  of  the  mixing  zone  for  Rayleigh-Taylor  unstable  interfaces  [16],  [  21].  Front 
tracking  has  provided  the  best  and  most  extensive  computations  for  this  problem  to  date.  In  contrast, 
modern  finite  difference  methods  without  interface  methods  have  yielded  values  about  30%  too  low  in  com¬ 
parison  to  experiment  for  the  growth  rate  [36].  In  spherical  geometry,  they  display  severe  mesh  orientation 
effects  as  well  [1].  Front  tracking  computations  were  used  in  systematic  studies  of  single  and  multi-mode 
Rayleigh-Taylor  interactions,  as  shown  for  example  in  Figure  1.  to  establish  the  dynamics  of  elementary 
multiphase  configurations  [21]  of  interacting  bubbles.  The  same  computations  are  now  being  used  to  test 
turbulence  models  and  multiphase  flow  models  in  the  high  Mach  number  compressible  context  [19]. 

Our  main  results  are  (a)  parallelization,  (b)  three  dimensional  tracked  computations  (in  progress),  (c) 
complex  wave  interactions  well  resolved  on  a  coarse  grid  and  (d)  the  use  of  interface  methods  in  tabular 
equations  of  state  for  multiphase  materials. 

A.  Parallelization.  The  parallelization  of  the  purely  hyperbolic  component  of  the  two-dimensional 
front  tracking  code  has  been  fully  implemented  on  the  INTEL  iPSC/^60  hvpercube.  enabling  the  parallel 
computation  of  gas  dynamics  problems.  This  parallelization  was  achieved  by  domain  decomposition  [13]. 
[  12].  [  11],  [  14].  The  spatial  domain  is  divided  into  a  union  of  disjoint  rectangular  subdomains,  with  the 
accompanying  division  of  the  tracked  physical  discontinuity  curves  among  the  subdomains.  An  extended 
boundary  region  of  n  mesh  blocks  in  each  direction  surrounds  each  subdomain  providing  overlap  into 
neighboring  subdomains.  Typically  n  is  an  upper  bound  for  the  finite  difference  method's  stencil  radius, 
for  example  n  =  2  for  the  Lax-Wendroff  method.  Thus,  the  boundary  region  for  the  i.j  mesh  block  lies 
entirely  in  the  eight  mesh  blocks  surrounding  it  (neglecting  the  slight  complication  of  physical  boundaries). 
In  Figure  l.  we  show  a  typical  interface  for  a  complex  fluid  mixing  process,  decomposed  into  16  subdomains, 
with  the  overlapping  boundaries  displayed  as  well. 

The  tracking  algorithm  progresses  iteratively  as: 
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Figure  1:  A  typical  interface  at  later  time  on  a  two-dimensional  320  X  300  mesh  is  decom¬ 
posed  into  16  sub-interfaces.  The  dotted  vertical  strips  shown  here  are  the  overlapping 
border  domains  shared  by  neighboring  processors. 


1.  Each  subdomain  updates  those  boundary  regions  of  neighboring  subdomains  that  lie  inside  its  own 
area.  (A  subdomain  never  updates  its  own  boundary  region,  which  lies  entirely  within  neighboring 
subdomains.) 

2.  The  discontinuities  are  propagated  and  solution  obtained  for  the  next  time  step  solely  within  each 
subdomain,  using  boundary  data  from  its  boundary  region,  which  is  stored  on  local  memory. 

This  scatter-gather  type  of  algorithm  is  most  efficient  if  all  subdomains  are  equally  busy  during  step  2. 
This  requires  a  sub-domain  assignment  algorithm  which  will  produce  subdomains,  not  of  equal  measure 
in  spatial  volume  but  of  equal  measure  in  propagation/solution  update  work.  At  present,  a  rather  simple 
sub-domain  assignment  is  in  place.  Current  gas  dynamics  calculations  have  shown  the  parallelized  code  to 
be  running  at  an  efficiency  of  approximately  90%.  The  scientific  power  can  be  seen  from  the  observation 
that  the  full  parameter  study  [8]  from  which  Figure  1  is  extracted  would  have  taken  an  estimated  17  years 
to  complete  on  a  Sun  Microsystems  SPARCstationl. 

B.  Three  dimensional  computations.  Of  critical  physical  interest  is  the  ability  to  do  computations 
in  three  spatial  dimensions.  We  have  begun  this  development  for  the  front  tracking  code.  Algorithms 
for  triangulating  surfaces  in  three  dimensions  and  for  re-gridding  dynamical  points  of  the  surfaces  have 
been  implemented.  A  preliminary  algorithm  for  generating  volume  filling  grids  which  match  moderately 
complicated  surfaces  was  constructed  but  further  work  in  this  area  is  required.  Additionally,  data  structures 
and  code  changes  to  handle  arbitrary  (spatial)  dimensionality  have  been  implemented  into  the  front  tracking 
code  to  support  calculations  in  1-,  2-  and  3-dimensions. 

C.  Complex  wave  interactions.  The  analysis  of  the  transitions  from  regular  to  irregular  refractions 
of  shock  waves  through  material  interfaces  has  resulted  in  an  improved  understanding  of  this  process.  The 
ideas  developed  here  were  applied  to  the  passage  of  a  shock  wave  through  a  bubble  [27],  yielding  a  substan¬ 
tial  improvement  in  numerical  resolution  of  the  refractions.  Figure  2  snows  illustrates  the  production  of 
an  anomalous  wave  refraction  during  the  shock-bubble  interaction.  As  a  result  of  this  work,  it  was  shown 
that  the  analysis  of  elementary  waves  given  in  [20]  was  incomplete  due  to  an  overly  restrictive  genericity 
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(a)  time  0.0  jisec 


Figure  2:  The  collision  of  a  shock  wave  in  water  with  an  air  bubble.  The  fluids  ahead  of 
the  shock  are  at  normal  conditions  of  1  atm.  pressure,  with  the  density  of  water  1  g/cc  and 
air  0.0012  g/cc.  The  pressure  behind  the  incident  shock  is  10  Kbar  with  a  shocked  water 
density  of  1.195  g/cc.  The  grid  is  60  X  60.  The  contours  in  Figure  2.b  and  Figure  2.c  are 

a  SCak  °f  °‘°01  ~  10  KbarS’  While  the  Pressure  ^nge  in  Figure  fg:anom.d  is  0 
ars.  The  tracked  fronts  are  shown  in  a  dark  line  superimposed  on  the  pressure 

contours. 


assumption,  and  that  at  least  one  additional  elementary  node,  the  total  internal  reflection,  must  be  al¬ 
lowed.  Subsequent  work  has  applied  these  ideas  to  the  Richtmyer- Meshkov  problem  of  a  shock  accelerated 

.  U/d  erface  f19]-  was  argUed  f26[  that  for  the  transition  of  a  regular  (self-similar)  shock  refraction 
m  o  an  irregular  configuration,  transitions  can  be  classified  into  five  basic  cases  depending  on  the  shock 
impedance  across  the  interface  and  on  whether  the  reflected  wave  is  a  shock  or  rarefaction.  This  classifi¬ 
cation  provides  the  background  necessary  to  incorporate  the  complex  configurations  produced  by  irregular 
shock  refractions  into  the  front  tracking  method.  Figure  3  illustrates  one  such  complex  interaction  that 
occurs  for  a  slow-fast  interaction  of  a  shock  with  a  fluid  interface.  Here  we  see  that  after  transition  the 
original  transmitted  wave  moves  ahead  of  the  incident  shock,  leading  to  a  complex  cascade  of  secondary 
wave  interactions.  Note  that  the  detail  of  wave  interaction  resolved  with  front  tracking  in  one  or  two  grid 
blocks  might  take  close  to  one  hundred  grid  blocks  for  comparable  resolution  by  other  methods.  The  anal¬ 
ysis  and  numerical  implementation  of  the  scattering  behavior  of  irregular  shock  refractions  is  an  essential 
component  for  full  scale  front  tracking  simulations  of  shock  accelerated  interfaces.  The  late  time  mixing 

e.  ?  -u6  mtTfaCe  1IUtlally  shown  in  Figure  3  is  §lven  ^  Figure  4.  It  should  be  noted  that  at  this 
point  the  fluid  interface  is  strongly  affected  by  shock  reflections  from  the  nearby  wall  at  the  bottom  of  the 
computational  domain. 

D.  Interface  Methods  for  Tabular  Equations  of  State.  The  front  tracking  software  was  used 
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Full  Computation 


Figure  3:  The  passage  of  a  shock  wave  through  a  random  interface  separating  two  gases 
of  differing  densities.  The  computation  is  shown  together  with  two  levels  of  graphical 
enlargement.  At  transition,  the  transmitted  shock  is  moving  faster  than  the  incident  shock 
leading  to  the  production  of  a  precursor  wave.  In  the  finest  enlargement,  one  can  see  complex 
wave  diffraction  patterns  resolved  down  to  the  level  of  a  single  mesh  block,  displaying  a 
unique  capability  of  front  tracking.  The  long  time  behavior  of  this  solution  (Figure  4) 
exhibits  interface  instability  similar  to  that  of  Figure  1. 


as  an  interpolation  scheme  for  piecewise  smooth  but  discontinuous  data.  This  method  was  applied  to 
the  representation  of  data  in  EOS  tables  with  phase  transitions  [10].  Starting  from  original  data  given 
on  a  coarse  grid,  spline  interpolation  was  used  as  an  initialization,  to  give  piecewise  smooth  data  defined 
on  a  finer  grid.  Lower  order  but  computationally  more  efficient  linear-bilinear  interpolation  was  used  to 
interpolate  the  functions  on  the  fine  grid.  This  mapping  of  data  from  coarse  to  fine  grid  is  computationally 
expensive,  but  is  only  performed  once.  For  repeated  evaluations,  a  significant  improvement  in  the  quality 
of  the  interpolated  data  was  obtained  in  this  way,  see  Figure  5.  This  interpolated,  piecewise  smooth  data 
was  then  used  to  solve  Riemann  problems  for  gasses  with  a  real  equation  of  state.  It  was  found  that  a  real 
gas  EOS  Riemann  problem  could  be  solved  in  no  more  than  about  3  to  8  times  the  time  required  for  a  7  - 
law  gas.  This  efficiency  depended  on  the  use  of  precomputed  and  preinverted  tables  for  the  sound  speed, 
the  Riemann  invariants  and  for  thermodynamic  variables  expressed  as  a  function  of  various  combinations 
of  independent  variables.  With  the  use  of  additional  tables,  a  time  at  most  3  times  the  7  law  gas  could 
have  been  achieved  in  all  cases  studied  [35]. 
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Figure  4:  A  detail  showing  the  late  time  chaotic  mixing  resulting  from  the  acceleration 
of  a  fluid  interface  by  a  shock  wave.  The  figure  shows  the  late  time  mixing  zone  for  the 
simulation  described  in  Figure  3.  The  detail  is  taken  from  the  bottom  of  the  computational 
cell  and  shows  the  shape  of  the  fluid  interface  from  Figure  3  500  microseconds  after  the 
shock  collision. 


3.2.2  Chaotic  Flows 

Our  results  provide  a  growing  body  of  knowledge  for  the  Rayleigh- Taylor  mixing  layer,  with  agreement 
among  theory,  experiment  and  direct  simulation  for  nearly  incompressible  flows.  We  also  found  a  surprising 
new  phenomenon,  discovered  computationally,  for  compressible  flows,  in  a  substantially  increased  growth 
rate  of  the  mixing  layer,  and  a  loss  of  universality  which  characterized  the  incompressible  case. 

The  mixing  process  is  characterized  by  a  penetration  height  h(t),  which  measures  the  distance  of  the 
most  advanced  bubble  from  the  position  of  the  unperturbed  interface.  This  height  obeys  a  scaling  law 
h  =  at2,  in  which  the  growth  constant  is  found  to  be  a  universal  constant  a  =  .06  in  the  incompressible 
experiments  of  Read  and  Youngs.  We  found  agreement  with  this  experimental  growth  rate  in  our  nearly 
incompressible  computations  [21]  and  [22],  but  found  surprising  new  phenomena  for  even  moderately 
compressible  fluids,  namely,  the  mixing  rate  a  may  exceed  twice  its  incompressible  value.  In  addition, 
some  loss  of  universality  was  observed  through  dependence  of  the  growth  rate  on  details  of  specification 
of  the  ensemble  of  initial  conditions  [8],  [  7].-  The  dependence  of  the  growth  rate  a  on  dimensionless 
compressibility  M2  is  shown  in  Figure  6. 

In  [23],  [  24],  and  [25]  a  statistical,  chaotic  theory  of  the  Rayleigh- Taylor  mixing  layer  is  given  in  terms  of 
a  renormalization  group  fixed  point  model.  The  renormalization  approach  is  used  since  the  chaotic  mixing 
layer  involves  dynamically  changing  length  scales.  The  model  is  validated  by  comparing  the  predicted 
growth  rate  of  the  mixing  layer  with  experiments  and  numerical  computations  [37].  The  three  main 
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(a)  (b) 


log  p  (g/cc)  log  p  (g/cc) 


Figure  5:  A  comparison  of  isobars  for  tabulated  equations  of  state.  The  plot  on  the  left 
shows  the  computed  isobars  for  a  SESAME  [31]  equation  of  state  near  a  phase  boundary. 
Since  the  rational  interpolator  used  in  SESAME  does  not  account  for  the  discontinuous 
derivatives  near  the  phase  boundary,  the  resulting  isobars  are  inaccurately  represented  in 
this  region.  The  plot  on  the  right  shows  the  computed  isobars  obtained  by  the  interface 
based  interpolator  which  accurately  represent  the  equation  of  state  properties  of  the  fluid 
near  phase  transition.  Both  plots  are  shown  with  the  phase  boundary  curve  superimposed 
as  a  dotted  line  over  the  isobars. 


ingredients  of  this  fixed  point  analysis  are  (a)  a  superposition  hypothesis  to  specify  the  bubble  -  bubble 
interaction  dynamics,  (b)  a  theory  of  single  bubble  dynamics,  (c)  a  statistical  model  to  incorporate  the 
above  solutions  to  the  one  and  two  body  problems  for  the  bubble  dynamics. 

The  superposition  hypothesis  was  given  in  [21]  and  [22].  It  describes  the  motion  of  the  outer  envelope 
in  the  Rayleigh-Taylor  mixing  layer  as  a  superposition  of  individual  bubble  envelopes.  The  superposition 
theory  explains  the  phenomenon  that  the  velocity  of  a  more  advanced  bubble  in  a  multi-bubble  system 
exceeds  the  velocity  of  a  bubble  of  the  same  size  in  the  single  bubble  system.  This  theory  also  explains  the 
fact  that  a  less  advanced  bubble  in  a  multi-bubble  system  changes  its  direction  of  movement  at  the  end  of 
its  period  of  interaction  with  a  more  advanced  neighboring  bubble.  The  superposition  hypothesis  has  the 
virtue  of  containing  no  free  parameter.  The  predictions  of  the  superposition  hypothesis  are  confirmed,  to 
within  the  accuracy  of  the  experiments,  for  the  incompressible  case  by  analysis  of  the  experiments  of  Read. 
Comparison  with  numerical  computations  for  compressible  fluids  shows  agreement  with  the  superposition 
hypothesis  for  small  compressibility  values,  but  reveals  disagreement' for  larger  values  of  compressibility, 
outside  the  range  covered  by  experiments.  An  explanation  for  the  disagreement  and  a  possible  basis 
for  modification  of  the  superposition  hypothesis  is  given  in  terms  of  density  stratification  of  the  fluids. 
Disagreement  is  also  noted  in  the  cases  where  bubble  splitting  occurs,  presumably  due  to  omission  of  high 
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Figure  6:  Plot  of  the  growth  rate  constant  vs.  the  compressibility.  The  vertical  bars  indicate 
the  variance  associated  with  the  choice  of  random  number  seed. 


frequency  bubble  splitting  modes  in  the  envelope  description.  Unexplained  disagreement  for  computations 
at  small  Atwood  number  is  also  noted. 

The  superposition  theory  depends  critically  on  a  good  description  of  the  single  (periodic)  bubble  dy¬ 
namics.  An  extension  of  a  previous  theory  for  the  growth  of  a  single  bubble  with  periodic  boundary 
conditions,  from  a  three  parameter  ODE  to  a  four  parameter  ODE.  was  presented  to  remove  an  earlier 
ansatz  which  lacked  physical  basis.  Two  of  the  four  parameters  of  the  new  theory  are  determined.  The 
remaining  two  must  be  determined  through  explicit  numerical  calculations.  Such  a  determination,  over  a 
limited  range  of  the  independent  variables  (Atwood  number  A  and  dimensionless  compressibility  A/2)  was 
presented  [22]. 

The  statistical  model  on  which  the  renormalization  group  fixed  point  is  based  describes  an  ensemble 
of  bubbles  of  the  same  radius,  whose  heights  are  defined  by  a  uniform  probability  measure  restricted  to  a 
bounded  interval.  The  statistical  dynamics  of  flow  with  bubble  merger  is  developed  by  treating  pairwise 
interactions  by  drawing  two  adjacent  bubbles  randomly  from  the  ensemble.  The  dynamics  of  each  pairwise 
merger  is  given  by  the  superposition  hypothesis  of  [21],  namely  that,  before  merger,  each  bubble  moves  with 
a  velocity  given  as  the  sum  of  a  scaled  single-bubble  velocity,  as  treated  in  [22].  and  an  envelope  velocity. 
The  bubble  of  higher  height  doubles  in  size  and  the  lower  bubble  is  removed  from  the  statistical  ensemble 
at  the  end  of  the  merger.  Differential  equations  are  then  obtained  for  the  common  radius,  average  height 
and  variance  of  height  of  the  ensemble  of  bubbles  as  a  function  of  time.  The  variance  of  height  is  shown 
to  have  a  natural  interval.  Its  lower  limit  is  a  trivial  fixed  point  corresponding  to  an  (unstable)  interface 
consisting  of  bubbles  of  identical  height.  Its  upper  endpoint  is  defined  by  instantaneous  merger  for  bubble 
pairs  of  extreme  separation.  By  studying  the  behavior  of  the  rate  of  change  of  variance  with  time  at  these 
two  endpoints,  the  existence  of  a  non-trivial  fixed  point  is  shown.  Figure  7  shows  a  numerical  verification 
of  this  renormalization  group  fixed  point  for  the  Rayleigh-Tavlor  instability. 

An  extensive  body  of  experiment  and  computation  predicts  a  constant  acceleration  for  the  leading  edge 
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Figure  7:  The  numerical  verification  of  a  renormalization  group  fixed  point  for  the  bubble 
envelope  of  a  Rayleigh-Taylor  unstable  interface.  The  graphs  shown  here  represent  the 
superposition  of  distinct  time  steps.  Both  axes  are  scaled  by  the  renormalization  group 
dynamics.  The  fixed  point  of  the  bubble  envelope  is  shown  in  each  graph.  For  the  second 
correlation,  on  the  right,  the  entire  graph  is  approximately  fixed  in  scaled  variables. 


of  this  mixing  region,  consistent  with  the  conclusions  of  the  fixed  point  predicted  by  this  theory.  The 
upper  and  lower  limits  placed  on  the  value  of  the  fixed  point  in  this  theory,  are  shown  to  yield  upper  and 
lower  limits  for  this  constant  acceleration  that  are  in  full  agreement  with  experiments  and  computations 
on  incompressible  and  nearly  incompressible  systems.  Further  studies  of  the  theory,  including  prediction 
of  transient  behavior,  dependence  on  density  ratio  and  compressibility,  assumptions  on  uniform  bubble 
radius,  and  extension  to  three  dimensions  remain  to  be  carried  out. 

We  are  now  studying  the  interior  of  the  mixing  zone  itself,  using  computational  data  from  well  resolved 
direct  simulation.  Statistical  analysis  of  fluctuating  quantities  reveals  struct  lire  which  is  more  complex  t  han 
simple  diffusion  [8],  [  7].  In  particular,  steady  acceleration  (Rayleigh-Taylor  unstable)  induced  mixing  of 
a  randomly  perturbed  interface  shows  non  monotone  density  contours  and  interior  structure  in  second 
order  correlations  of  fluctuating  quantities.  This  is  consistent  with  theories  of  turbulent  boundary  layers, 
which  show  at  least  three  distinct  regions  within  the  mixing  layer.  The  study  of  fluctuating  quantities 
suggests  renormalization  group  fixed  point  behavior  near  the  edge  of  the  mixing  zone  only  and  in  the 
nearly  incompressible  case  only  [8],  [  7], 

3.2.3  Nonlinear  Waves  and  Nonlinear  Materials 

Striking  new  developments  in  the  theory  of  nonlinear  waves  for  conservation  laws  have  given  rise  to  a  now 
picture  for  wave  interactions  in  the  large.  Improved  computational  algorithms  and  modeling  of  physical 
phenomena  are  to  be  expected  as  the  consequences  of  these  results  are  explored.  It  is  now  widely  recognized 
that  nonlinear  waves  may  contain  two  or  more  significant  length  scales.  As  a  consequence  the  ratios  of 
these  length  scales  become  important  dimensionless  parameters,  controlling  macroscopic  wave  speed  and 
structure  [29]  for  three  phase  flow,  chemically  reactive  flow,  elasticity  and  MHD  waves  [3], 

Plohr.  Marchesin  et.  al.  [28]  have  produced  a  very  important  unifying  framework  for  the  fundamental 


Figure  8:  Shown  here  is  a  two  dimensional  slice  of  the  wave  manifold  W  for  a  representative 
2x2  system  of  conservation  laws  with  an  elliptic  region.  The  figure  shows  the  influence  of 
Hopf  bifurcation  on  the  admissibilty  of  shock  waves,  when  the  viscous  admissibility  criterion 
is  taken  into  account.  The  connecting  orbit,  oriented  from  U-  to  U+,  represents  a  traveling 
wave  between  these  two  states,  i.e.  a  viscous-admissible  shock  wave.  The  figure  shows  a 
“bubble”  region  in  which  the  connection  between  U-  and  U+  is  prevented  by  the  occurence 
of  a  limit  cycle,  which  starts  at  the  Hopf  manifold.  The  limit  cycle  ends  in  the  homoclinic 
loop,  beyond  which  the  connection  from  U _  to  U+  is  established.  The  live  phase  diagrams, 
with  associated  singular  points  {/_  and  t/"+,  show  orbits  in  the  three  open  regions  to  the 
right  of  the  characteristic  manifold,  and  at  two  transitional  boundary  points.  Only  the 
shocks  in  the  shaded  region  are  admissible. 


waves  occurring  in  general  systems  of  n  conservation  laws.  By  using  a  local  change  of  coordinates,  the 
Rankine-Hugoniot  relation  is  shown  to  take  the  form  1Z  ■  T  =  0.  The  solution  set  where  1Z  =  0  is  a 
trivial  solution  set  representing  constant  states  and  rarefaction  waves.  By  eliminating  this  trivial  solution 
set,  it  is  shown  that  the  solution  set  >V,  where  T  =  0,  is  a  smooth  manifold  of  dimension  n  +  1  and  is 
the  closure  of  the  set  of  shock  points.  W  is  termed  the  fundamental  wave  manifold.  Significantly,  both 
rarefaction  and  shock  waves  are  represented  within  >V,  in  accordance  with  the  heuristic  idea  that  shocks 
of  infinitesimal  strength  are  infinitesimal  rarefaction  fans.  The  rarefaction  points  form  an  re-dimensional 
submanifold  C  of  W,  the  characteristic  manifold,  which  may  have  singularities  at  points  corresponding  to 
coincident  wave  speeds.  The  familiar  rarefaction  curves  in  state  space  for  the  system  of  conservation  laws  are 
projections  of  a  single  family  of  curves  in  C  which  form  a  one-dimensional  foliation  of  C.  Correspondingly, 
the  manifold  W  is  foliated  by  two  families  of  curves,  called  shock  curves  since  they  project  onto  the  classical 
shock  curves  in  state  space.  The  work  shows  how  this  wave  manifold  framework  can  shed  light  on  two 
fundamental  problems  in  the  theory  of  conservation  laws.  The  first  is  the  physical  admissibility  of  shock 
waves  determined  by  properties  of  dynamical  systems  parameterized  by  the  points  of  VV.  The  second 
problem  is  the  bifurcation  of  wave  curves,  which  correspond  to  loss  of  transversality  between  rarefaction, 
shock,  and  composite  foliations,  and  to  the  boundary  of  the  region  of  admissible  waves. 

The  wave  manifold,  VV,  contains  many  nonadmissible  shock  waves.  The  most  fundamental  notion  of 
admissibility  presently  known  is  the  viscous  profile  criterion,  which  states  the  shock  wave  must  be  the  limit, 
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as  the  viscosity  tends  to  zero,  of  traveling  waves  for  an  associated  viscous  conservation  law.  This  viscous 
conservation  law  gives  rise  to  a  dynamical  system  with  critical  points  corresponding  to  the  states  on  the 
left  and  right  of  the  shock  wave.  In  [6],  some  of  the  mathematical  issues  which  will  be  associated  with  loss 
of  admissibility  are  studied.  Of  particular  interest  is  the  demonstration  that  Hopf  bifurcation  can  be  the 
mechanism  which  leads  to  loss  of  admissibility.  Similarly,  homoclinic  orbits  are  associated  with  a  loss  of 
admissibility.  In  Figure  8  we  see  that  the  connecting  orbit  bifurcates  when  crossing  the  Hopf  bifurcation 
locus,  so  that  one  end  is  connected  to  the  limit  cycle  emerging  from  the  Hopf  bifurcation.  The  connection 
between  the  end  states  of  the  shock  wave  is  thus  broken,  implying  the  nonadmissibility  of  this  shock  wave. 

The  transitional  waves  are  the  most  curious  of  the  novel  shock  waves  discovered  in  the  recent  renewal  of 
interest  in  Riemann  problems.  These  waves  have  dynamical  system  orbits  which  connect  saddle  points  to 
saddle  points,  and  thus  they  appear  to  be  inherently  unstable.  However,  they  have  been  essential  to  obtain 
a  satisfactory  existence  and  uniqueness  theory  for  solutions  of  Riemann  problems.  Stability  analysis  has 
been  used  in  the  search  for  a  more  satisfactory  basis  for  accepting  these  waves  as  physically  meaningful, 
In  [38]  nonlinear  stability  was  studied,  and  on  a  numerical  level,  established  for  these  shock  waves. 

Modeling  of  phase  transitions  was  shown  to  depend  on  an  additional  degree  of  freedom  (the  order 
parameter),  enlarging  the  system  to  give  internal  structure  [17].  T.-P.  Liu  in  his  earlier  work  on  nonequi¬ 
librium  thermodynamics  [30]  also  considered  a  larger  system  as  a  regularization  of  a  smaller  one.  Related 
studies  of  the  phase  field  model  of  phase  boundaries  [4],  [  5]  follow  this  point  of  view.  The  relation  between 
kinks  and  the  loss  of  convexity  in  the  equation  of  state  and  multiple  or  split  waves  in  Riemann  problems 
was  developed  in  a  systematic  fashion  in  [32]. 

A  fully  conservative  Eulerian  formulation  of  the  elasticity  equations  [33],  [  34]  has  recently  been  ob¬ 
tained.  This  formulation  promises  to  be  of  considerable  importance.  In  cases  of  large  deformation,  Eulerian 
computations  are  necessary  to  avoid  the  severe  mesh  distortion  of  Lagrangian  grids.  A  thermodynamically 
consistent  form  of  the  elasticity  constituitive  laws  derived  from  a  free  energy  with  a  small  deviatoric  strain 
was  given  [33],  [  34].  This  formulation  allows  arbitrary  fluid  behavior  in  the  pressure  and  thermal  modes, 
and  depends  on  a  single  shear  modulus  to  describe  shear  strength.  It  was  shown  that  this  formulation 
is  the  lowest  order  approximation  to  a  general  free  energy  in  the  case  of  a  small  deviatoric  strain.  An 
elasticity  Riemann  solver  was  constructed  for  this  case.  Computations  have  shown  the  distinct  advantage 
of  the  fully  conservative  approach  [2],  [  15]. 

In  our  study  of  nonlinear  materials,  a  new  and  fully  conservative  formulation  of  plasticity  was  also 
formulated  [34].  Based  on  the  importance  of  the  conservation  formulation  for  computations  in  gas  dynamics, 
we  expect  that  this  discovery  could  be  of  very  fundamental  character. 
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